# compare with ~/xchg/Match/cem/a1lalonde_nocaliper.R, except with estimand=ATT

#the truth for lalonde, regression model

#compare with lapo:~/xchg/Match/GenMatch/mc1/lalonde3matching1.Rout

library(Matching)
data(lalonde)

sims <- 1000
estimand <- "ATT"
obs <- nrow(lalonde)

#setDATA seed
set.seed(2810192)

attach(lalonde)
propensity  <- glm(treat~ I(age^2) + I(educ^2) + black +
                   hisp + married + nodegr + I(re74^2) + I(re75^2) +
                   u74 + u75, family=binomial, data=lalonde)

X <- cbind(rep(1, length(treat)),
           propensity$linear.pred,
           I(log(age)^2),
           I(log(educ)^2),
           I(log(re74+0.01)^2),
           I(log(re75+0.01)^2))

propensity.coeffs <- propensity$coeff
propensity.coeffs <- as.matrix(c(
                                 1+00,  #(Intercept)        
                                 .5,    #linear.pred
                                 .01, #age
                                 -.3,  #educ
                                 -0.01, #I(re74^2)          
                                 0.01   #I(re75^2)          
                                 ))

mu = X %*% propensity.coeffs
Tr.pred <- exp(mu)/(1+exp(mu))
#print(summary(Tr.pred))

TreatmentEffect <- 1000
TreatmentReal <- matrix(nrow=obs, ncol=1)

#return objects
raw <- matrix(nrow=sims)
correct.est.bias <- matrix(nrow=sims)
correct.est.nobias <- matrix(nrow=sims)
est.bias <- matrix(nrow=sims)
est.nobias <- matrix(nrow=sims)
iv.est.bias <- matrix(nrow=sims)
iv.est.nobias <- matrix(nrow=sims)
ps.est.bias <- matrix(nrow=sims)
ps.est.nobias <- matrix(nrow=sims)
psmh.est.bias <- matrix(nrow=sims)
psmh.est.nobias <- matrix(nrow=sims)


results <- function(s)
  {
    cat("RAW bias:", mean(raw[1:s]), "\n")
    cat("RAW RMSE:", sqrt(mean(raw[1:s]^2)), "\n")    
    
    cat("\n")
  }

sims.data <- vector("list",sims)
count <- 0
for (s in 1:sims)
  {

    cat("\nS:", s, "\n")
    count <-  count+1
    
    for(i in 1:obs)
      {
        TreatmentReal[i] = sample(0:1, 1, prob=c(1-Tr.pred[i],Tr.pred[i]))
      }

    Y <- I(TreatmentEffect*TreatmentReal) + .1*exp(.7*log(re74+0.01) + .7*log(re75+0.01)) + rnorm(obs, 0, 10)

    sims.data[[s]] <- cbind(Y, TreatmentReal)
    

    raw[s] = mean(Y[TreatmentReal==1])-mean(Y[TreatmentReal==0]) - TreatmentEffect
    cat("raw[",s,"]:", raw[s], "\n")
    cat("raw old[",s,"]:", (mean(Y[TreatmentReal==1])-mean(Y[TreatmentReal==0])) , "\n")
    
    if (count==10)
      {
        results(s)
        count <-  0
#        save.image(file=".RDATA.lalonde3matching1")
      }
  } #end of sims

results(s)

save(sims.data,file="RData.mc2.data1")
